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Abstract. We obtain empirical formulae for the final remnant black hole mass, spin, 

and recoil velocity from merging black-hole binaries with arbitrary mass ratios and 
spins. Our formulae are based on the mass ratio and spin dependence of the post- 
Newtonian expressions for the instantaneous radiated energy, linear momentum, and 
angular momentum, as well as the ISCO binding energy and angular momentum. The 
relative weight between the different terms is fixed by amplitude parameters chosen 
through a least-squares fit of recently available fully nonlinear numerical simulations. 
These formulae can be used for statistical studies of N-body simulations of galaxy 
cores and clusters, and the cosmological growth of supcrmassivc black holes. As an 
example, we use these formulae to obtain a universal spin magnitude distribution of 
merged black holes and recoil velocity distributions for dry and hot /cold wet mergers. 
We also revisit the long term orbital precession and resonances and discuss how they 
affect spin distributions before the merging regime. 
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1. Introduction 

Black holes at the centers of galaxies and globular clusters significantly impact the 
dynamical evolution of these astronomical structures. Of particular importance to the 
dynamics are the black-hole (BH) mass, spin, and location (if off-center); properties that 
can significantly change following a BH merger. When two galaxies merge, an event that 
is expected to occur several times during a galaxy's evolution, the supermassive BH at 
their centers form a black-hole binary (BHB) that eventually inspirals and merges. 
Similarly, intermediate-mass BHs in globular clusters can form tight BHBs that inspiral 
and merge. Consequently, accurate models for the mass, spin, and gravitational recoil of 
the merger remnants of BHBs are of great astrophysical interest. However, these models 
require accurate simulations of merging BHs, a problem in the highly-nonlinear regime, 
that only recently became feasible due to breakthroughs in Numerical Relativity [HEllS]. 

The first attempts at modeling the remnant BH of BHB mergers using fully 
nonlinear simulations utilized the 'Lazarus method' j4j, which combined short-term 
numerical simulations of BHBs, just prior to merger, with perturbative calculations. 
With the advent of the 'moving punctures' [21 E] and generalized harmonic [T] 
approaches, it became possible to accurately model merging BHBs from inspiral through 
merger and ringdown using fully-nonlinear numerical simulations. As a result of these 
breakthroughs, NR groups from around the world have been able to develop heuristic 
models for the properties of remnant BHs as a function of the orbital and intrinsic BH 
parameters of the binary (at-least for part of the parameter space). 

The initial attempts at modeling the properties of remnant BH focused on the 
mass and spin using ad hoc interpolation formulae. In [H], El Ej we studied equal-mass, 
spinning BHBs, where the individual BH spins were aligned and counter-aligned with 
the orbital angular momentum, using fully nonlinear numerical calculations. We found 
a simple quadratic polynomial relating the final mass and spin of the remnant with the 
spins of the individual BHs. This scenario was later revisited in [HI |9], and in [10] the 
authors generalized the formula for the remnant spin (by assuming that the angular 
momentum is only radiated along the orbital axis, and neglecting the energy loss) in 
order to model arbitrary BH configurations (these assumptions were relaxed in a follow- 
up paper [H]). A generic formulae for the final spin was proposed in [12] based on 
simulations with aligned and non-aligned spins. A more comprehensive approach, using 
a generic Taylor expansion consistent with the the physical symmetries of the problem, 
and with parameters chosen by a least-squares fit to many simulations, was developed 
in [13] . All of these models used low-order polynomial interpolation functions to predict 
the remnant mass and spin as a function of the individual BH masses and spins. On 
the other hand, in [T^, a different approach, based on approximate analytic models 
for the merger, was used. Here the authors extended the particle limit approximation 
for the radiated mass and angular momentum to the comparable- mass regime; ignoring 
effects of post-ISCO (Innermost Stable Circular Orbit) gravitational radiation. This 
approach was further improved in |15J by taking binding energies into account. All of 
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these approaches show a certain degree of agreement with the remnant masses and spins 
obtained in the few dozen fully nonlinear numerical simulations available, but significant 
uncertainties concerning accuracy outside this range of the parameter space remain. 
Here we propose a set of formulae that incorporate the benefits of both approaches and 
regimes in a unified way; using analytic techniques to develop empirical models with 
free parameters determined by numerical results. 

Due to its significant impact on astrophysics, the modeling of the remnant recoil 
followed an independent path, particularly since the discovery [16l HT] that the spins 
of the black holes play a crucial role in producing recoils of up to 4000 kms~^. 
The realization that the merger of BHBs can produce recoil velocities that allow the 
remnant to escape from major galaxies led to numerous theoretical and observational 
efforts to find traces of this phenomenon. Several studies made predictions of specific 
observational features of recoiling supermassive black holes in the cores of galaxies in 
the electromagnetic spectrum [HI [TH |20l |2ll [221 [231 |2l]. Notably, there began to 
appear observations indicating the possibility of detection of such effects [251 ESI [27], 
and although alternative explanations are possible [281 ISll |30] , there is still the exciting 
possibility that these observations can lead to the first confirmation of a prediction of 
General Relativity in the highly-dynamical, strong-field regime. 

This paper is organized as follows. In Sec. |2] we describe our empirical formula for 
the remnant gravitational recoil and provide the leading coefficients for this formula. In 
Sec. |3] we describe our formula for the final remnant mass, while in Sec. |4] we describe 
the formula for the final remnant spins. We provide fits to the constants in the remnant 
mass and spin formula in Sec. |5| In Sec. |6] we revisit the gravitational alignment and 
antialignment mechanisms for long term inspiral orbits, and discuss the consequences 
and apphcations of our formulae in Sec. [7| 



2. Remnant Recoil Velocities 



In our approach to the recoil problem [161 IE] we used post-Newtonian (PN) theory as 
a guide to model the recoil dependence on the physical parameters of the progenitor 
BHB (See Eqs. (3.31) in [3l|), while arguing that only full numerical simulations can 
produce the correct amplitude of the effect (see Eq. (|3| below). For example, in the 
instantaneous radiated linear momentum, there are terms of the form 



dP T] 



2 



■ ■ ■ + 



F{r,v) ■ {d2 - qdi] 



(1) 



dt l + q 

where L is the unit vector pointing along the instantaneous orbital angular momentum, 
F{f, v) is a vector in the orbital plane that is only a function of the orbital position and 
its time derivative, q = mi/m2 < 1 is the mass ratio, rj = q/{l + g)^ is the symmetric 
mass ratio, and cJj = Si/mf is the intrinsic spin on black hole i. We incorporate this 
term by adding a term 

Kecoil = ■ ■ ■ + (^(I^ 0«2 - QC^il COs(eA - 60)]^ t (2) 
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to our fitting formula for the recoil velocity (see Eq. ([s]) below), where the fitting 
constants K and Gq approximate the net effect of the dynamics of this term during 
the last few M of the rapid plunge (where most of the recoil is generated) and Ga 
is the angle that A = M^(a2 — '?«i)/(l + q) where M = mi + m2, makes with the 
infall direction at merger. Our heuristic formula describing the recoil velocity of BHB 
remnants was theoretically verified in several ways. In [17] we confirmed the sinusoidal 
dependence [cos Ga in Eq. ([s])] of the recoil on the direction of the in-plane spin for 
the so-called 'superkick' configurations, a result that was reproduced in [32] for binaries 
with different initial separations. While in [33] the authors verified the decomposition 
of the spin-dependence of the recoil into spin components perpendicular and parallel to 
the orbital plane. Similarly, in [31] the authors determined that the quadratic-in-spin 
corrections to the in-plane recoil velocity are less than 5% of the total recoil. Recently 
in [35] we confirmed the leading r]"^ (where r) is symmetric mass ratio) dependence of the 
large recoils out of the orbital plane (see also [36]). 

Here we augment our original empirical formula with subleading terms, higher order 
in the mass ratio, and include a new term linear in the total spin, motivated by higher 
order post-Newtonian computations 



Vrecoii{q,a) =Vmei + f_L(cos^ei + sin ^62) + v\\ n\\ 



v± = H 



1 + g) 
2 



K ^ 



1 + g) 
2 



i + q) 



(1 + Bh v) (4 - qal) + Hs (a| + q'al) 

1 + B^f]) I '^2' ~ q'^il cos(Ga — Go) 



'1 - q) 

+ Ks ,\ , ^ I + q^a^ | cos(G5 - Gi 



(3) 



[i + q) 

where the index ± and || refer to perpendicular and parallel to the orbital angular 

momentum respectively and n\\ = L. 61,62 are orthogonal unit vectors in the orbital 
plane, and ^ measures the angle between the unequal mass and spin contribution to the 
recoil velocity in the orbital plane. The new constants Hs and Ks can be determined 
from new generic BHB simulations as the data become available. The angles, Ga and 
G5, are the angles between the in-plane component of A = M{S2/m2 — Si/mi) or 
S = 81 + 82 and the infall direction at merger. Phases Gq and Gi depend on the 
initial separation of the holes for quasicircular orbits. A crucial observation is that 
the dominant contribution to the recoil is generated near the time of formation of the 
common horizon of the merging black holes (See, for instance Fig. 6 in [38j). The 
formula ^ above describing the recoil applies at this moment (or averaged coefficients 
around this maximum generation of recoil) , and has proven to represent the distribution 
of velocities with sufficient accuracy for astrophysical applications. 

The most recent estimates for the above parameters can be found in [35j and 
references therein. The current best estimates are: A = 1.2 X 10^ km s'\ B = -0.93, 
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H = (6.9 ±0.5) X 10^ km s~\ K = (6.0 ±0.1) x 10^ km s"\ and ^ ~ 145°. Note that we 
can use the data from [35] to obtain K = (6.072 ±0.065) x 10^ kms~^, if we assume that 
Bk and Ks are neghgible. Finally, if we fit the data to find K and Ks simultaneously 
we obtain K = (6.20 ±0.12) x 10^ kms"^ and Ks = -0.056 ±0.041, where we made the 
additional assumption that Gq = Oi (since S = A for these runs). An attempt to fit K, 
Ks, Bk simultaneously does not produce robust results with currently available data 
(one of the reasons for this is that different values of K and Bx produce very similar 
predicted recoil velocities over the range 0.1 < q < 1). Note that the values for the 
dominant K term are reasonably insensitive to the different choices for the fits, while 
finding the subleading terms require additional runs and higher accuracy. 

The above equation ([s]) contains all the expected linear terms in the spin, and 
includes ten fitting parameters. Based on the works [HZ] one could add quadratic terms, 
and this will be published elsewhere by the authors. 

From a practical point of view, for statistical simulations of BHB mergers, where 
the directions of the spins and infall direction is not known, one should take a uniform 
distribution for the in plane-components of di and 6:2 over all possible angles, define 
G5 and 0A with respect to a fixed arbitrary in-plane vector (say x), and take Oq = 0. 
The resulting distribution of recoil velocities will be independent of the choice of the 
arbitrary in-plane vector (but will depend weakly on 0i). If we ignore the subleading Ks 
correction, then 0i will not enter the recoil calculation. It's effects can be incorporated 
by including the Ks term and averaging over all possible values of Gi. 



3. Remnant Mass 



Motivated by the success of the empirical formula for the recoil, we propose a new 
empirical formula for the total radiated energy based on the post-Newtonian equations 
for the instantaneous radiated energy (see Eqs. (3.25) in [31], and for the quadratic 
terms in the spin see Eq. (5.4) in |37j): 

6M/M = r] Eisco + ^2^?' ± ^3^' 

"^ (1 g)2 ("2 + ± (1 - q) (a| -qa\) + Ea |«2 ± g^iP 

+Eb + ga^l^ (cos^(e+ - 62) ± Ec) + Ed |«2 - qai\^ 

+Ee - qa^l^ {cos^Q. - 83) ± ^f) |, (4) 

where G± are the angles that A± = M{S2/m2 ± Si/rrii) make with the radial direction 
during the final plunge and merger (for comparable-mass BHs, a sizable fraction of the 
radiation is emitted during this final plunge, see for instance Fig. 6 in [38]). Phases 02,3 
are parameters that give the angle of maximum radiation for these terms, and depend 
on the initial separation and parameters of the binary at the beginning of the numerical 
simulation. 

In addition to the terms arising from the instantaneous radiated energy, which gives 
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12 fitting parameters, we also included terms associated with the secular loss of energy 
in the inspiral period from essentially infinite separation down to the plunge. In order 
to model this contribution we adopted the effective one body form [39] supplemented by 
the rj"^ effects from self force calculations [IQ] and 2PN effects of the spins (see Eq. (4.6) 
in [31]), to obtain 

Eisco ^ (1 - V8/3) + 0. 1038037/ 

1 



36v^(l + g)2 
5 



q{l + 2q)al + {2 + q)al 
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\\J\ 



[a. 



l^2^ 



(5) 



324^2(1 + qy 
+0{a^). 

The above expression only includes quadratic-in-spin terms for compactness, hence it is 
expected to produce reliable results for intrinsic spin magnitudes < 0.8 (because 
the binding energy is a very steep function of a for a > 0.8 and the quadratic 
expressions above are no longer appropriate). Note that we used the full expressions 
from [31] to obtain our fitting parameters. Here we fit the leading-order parameters 
using available data, and as new data become available, we expect to be able fit the 
remaining parameters. 



4. Remnant Spin 



In an analogous way, we propose an empirical formula for the final remnant spin based 
on the post-Newtonian equations for the radiated angular momentum and the angular 
momentum of a circular binary at close separations (see Eqs. (3.28) and (4.7) in [21]), 

-2 



ttfinal 



6M/My 



7] 



(1 
+ (1 



qy 



r]Jisco + {J2V^ + J^n^) n\\ 
JA{a\ + q^a\) + Jb{1- 



q) {al -qa\) 



q) - ga^lv^JA cos[2(0A - 64)] + J\mn^ 
+ 14 + g'a^l VJ5Cos[2(05 - 65)] + Jusfi^ }. (6) 

Note that, even at linear order, there are important contributions of generic spinning 
black holes producing radiation in directions off the orbital axis that do not vanish in the 
equal-mass or zero-total-spin cases. The above formula can be augmented by quadratic- 
in-the-spins terms [S7l HI] of a form similar to the terms added to the radiated energy 
formula Q. However, these terms are less significant for modeling the final spin (see, 
for instance Fig. 21 of [7]). 

Again, we use the effective one body resummation form [39], supplemented with the 
rf effects from self force calculations [iQ] and the 2PN effects of the spins (see Eq. (4.7) 
in [21]), to obtain 

Jisco ~ I 2 - 1.52558627] - ^ 



9^2(1 + ^)2 



g(7 + 8g)ai + (8 + 7g)a^ 
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9^2(1 + qY'"'^ ' ' ' ' 7] {l + qY — V ; 

This expression represents a quadratic expansion in the spin-dependence, hence we 
expect to produce rehable results for intrinsic spin magnitudes a^ < 0.8 (hence 
ttfinal < 0.9). 



5. Determination of fitting parameters 

Here we show how resuhs from current full numerical simulations can be used to 
determine the fitting constants in the equations for the final remnant mass and spins 
of a BHB merger. This procedure can be repeated and extended as we have access 
to new runs and can also help in designing new simulations to optimally determine all 
fitting constants and better cover the 7-dimensional physical parameter space of BHBs. 
We used Mathematica's LinearRegression and NonLinearRegress functions to find the 
fitting parameters and estimate the errors in the parameters. Our method for finding 
the fitting parameters was to first fit to simulations with symmetries that caused most 
terms to vanish in order to fit to as few parameters at a time as possible. Then, after 
fixing the parameters we found in earlier fits, we fit to simulations with less symmetry 
to obtain other parameters. For example, we first find E2 and i?3, and then using these 
values, fit additional data to obtain Es-, etc. 

Energy radiated: For the non-spinning case, we fit the data from 8 simulations 
found in |l2l |l3] (see also [Hj). Here we fit -ERad versus 77, where E-^^a is the total 
radiated energy for a given configuration minus the binding energy of the initial 
configuration (where the binding energy is negative). We calculate the binding energy 
using the 3PN accurate expressions given in [13]. A fit of the resulting data gives 
E2 = 0.341 ± 0.014 and E3 = 0.522 ± 0.062. In order to estimate Es, Ea, Ea, and 
Ed, We use the remnant masses from 13 simulations for spinning BHBs with spins 
aligned with the orbital angular momentum given in [121 HTj (see also [5j), and find 
Es = 0.673 ± 0.035, Ea = -0.36 ± 0.37, Ea = -0.014 ± 0.021, and Ed = 0.26 ± 0.44. 
These large uncertainties in the fitting parameters are due to the effect of correcting 
for the binding energy in these simulations. Finally, fits from the final remnant masses 
from 5 simulations [T7] yields Ee = 0.09594 ± 0.00045 and fits from 5 equal-mass 
configurations in [35] yield Eb = 0.045±0.010. An accurate fit to Ed is not possible with 
the configurations available in [35]. Note that our fits for Ea, Ea, and Ed are consistent 
with the parameters set to zero. This is due to the fact that the errors introduced in 
renormalizing the data are of the same order as the effects of these subleading terms. 

Angular momentum radiated: For the non-spinning case, we fit the data from 8 
simulations in ^ |13], and find J2 = -2.81 ± 0.11 and J3 = 1.69 ± 0.51. A fit to Ja 
and Jb from 13 simulations in [IHIIIT] yields Ja = -2.97±0.26 and Jb = -1.73 ±0.80. 
However, we determined that the uncertainty in J a and Jb is actually closer to 1.0 by 
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considering fits to the independent datasets in |46j. 

From the combined fit, we find that 2.42% < 5M/M < 9.45% and 0.34 < agnai < 
0.92 for the equal-mass, aligned spin scenario, in the region where the fit is vahd 

(lafinall < 0.9). 

Finally, we note that much of the errors in the fitting parameters are due to 
differences in the normalizations between the various runs. Some authors choose 
normalize their simulations such that mi + m2 = 1, which approximates a binary that 
inspiraled from infinity with an initial mass of 1, while others choose to normalize their 
simulations such that the initial ADM mass is 1. In this latter case we attempted to 
renormalize the results using the 3PN expression for the binding energy. However, the 
errors introduced by renormalizing data, or assuming that the ADM mass at infinite 
separation is 1, introduces uncertainties in our fitting parameters. This affects both 
6M/M directly and cJfinai indirectly through 6M/M. Ideally we would use a set of 
simulations with the same normalization and all starting from the same initial orbital 
frequency. 

From a practical point of view, for statistical simulations of BHB mergers, where the 
infall direction and the directions of the spins at merger are not known, one should take 
a uniform distribution for the in plane-components of di and a2 over all possible angles, 
define the angles G5, Ga, and 0+ (note 0_ = ©a) with respect to a fixed arbitrary 
in-plane vector (say x), and take a uniform distribution for the unknown angles Oi.3,5. 
The angles 60,2,4 can be set to zero, since the final distributions will be independent of 
this choice (the distribution will only be a weak function of the relative angles Go — Oi, 
etc.). The resulting distributions will be independent of the choice of the arbitrary in- 
plane vector (but will depend weakly on 61^3^5). However, the angles ©1,3,5 only appear 
in subleading expressions and the uncertainties in the final distributions of the spins, 
masses, and recoils should not be significant for astrophysical applications. 

6. Inspiral phase 

One of the important application of our formulae is to study statistical distributions 
of the final mass, spin and recoil of the remnant merged black hole given an initial 
distribution of individual spins and mass ratios. This kind of studies have been 
performed lately, see for instance [IE] , assuming initial random distribution of individual 
spin directions and magnitudes as well as mass ratios. This choice was supported by 
the post-Newtonian studies |1S] that in the (dry) inspiral phase, preliminary to the final 
merger we have modeled in this paper, there is not an strong alignment of the spins 
with the orbital angular momentum, as there would be if, for instance, we would have 
large accretion of gas in the system (wet mergers). 

The simulations in [38] actually found that, gravitational radiation induced 
precession of the orbital plane during the inspiral phase leads to an small bias of the 
spins towards counter- alignment. These results were the product of integrating 3.5PN 
equations of motion form separations r = 50M down to the merger regime around 
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r = 5M. It was point out by studying averaged PN equations in the quasicircular 
orbits regime |^ [50| |5T] , that on longer time scales there are resonances that might 
affect the distributions of spin directions by the time of merger. Since these studies are 
complementary to those presented in |18], we will investigate this issue analytically at 
a lower PN order than in the numerical studies of Ref . [Hj , but retaining the radiation 
reaction effects on the orbital plane for consistency with the integration of the PN 
equations of motion in the Hamiltonian formalism. 

In terms of the notation and approach of Ref. |18] we consider 

, T* , L' S-j L ' S-j L ' Sj \L\ L ' Sj ISA 

(L ■ SA = -^-^ + ^ I ' , (8) 

irllQI irllQI irlsiQI |r||c|2 ^ ' 

where we can set 15*41' = at this order of approximation. In |1S] Si and the conservative 

part of L terms did not contribute due to the nature of the statistical studies performed 
in that paper. We hence focus on the dissipative part. 

(T B \ _ -^dis - Si L ■ Si |I/|dis , . 



With the PN techniques described |18] we find 

(f - ^ ^" g ^ 

^ ^ ^''"'^ 15 M {l + qY |ai| 

X |g (61g + 48)(P.ai)2 + (61 + 48g)(P.ai)(p.52)|, (10) 

- - . ^ v}} q 1 

{L ■ 5*2)^18 ^ 



15 M {l + qY lasi 

X |g (61g + 48)(P-ai)(P-a2) + (61 + 48g)(P-a2)'}. (11) 

Note that this expressions are defined negative when averaged over spin directions with 
only the squared terms (P ■ diY contributing. By integrating them over time, we obtain 

similar results to the expression for (L • Sy^:^^ in Eq (18) of |18], that lead us to the 
conclusion that distributions of spins show some bias towards counter-alignment with 
respect to the orbital angular momentum. Note that the instantaneous counteralignment 
mechanism acts at every radius, with increasing strength for small separations, where 
the orbital velocity is large. 

To investigate the small mass ratio limit, i.e. g — )■ 0, we compute the time integral 



(roughly speaking, multiply by 1/q) of Eq. (11), for instance. We can then see that if 
the larger black hole's spin, 5*2, is initially randomly distributed in the limit g — )■ it 
ends up with some counteralignment. On the other hand, the smaller black hole's spin, 
5*1, would remain to be random oriented as seen from the vanishing of the right hand 



side of Eq. (jlOj). 

Note that the above equations do not use orbital averages since the effect is 
particularly strong in the latest part of the inspiral, when averages are not a good 
approximation. In the alternative regime, when the inspiral motion is very slow. 
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resonance orbits have been found using orbit averaged descriptions [IHl Eni EH 152] . 
These resonance orbits lead to ahgnment or antiahgnment of spins if one starts from 
an initial aligned or antialigned large hole and allow random orientations for the less 
massive one. Note that if both spins are allowed to be chosen at random initially, as 
we assumed in our computations, then the resulting evolution leads to still random 
distributed spins. 

The resonance mechanism is complementary to the mechanism we studied in [15] . 
The former takes place on very long time scales compared to precession, while the later 
mechanism is quadratic in the spins (as seen in Eqs. (10) and (11 ), hence higher order.) 
In order to quantify which of them is the predominant mechanism long term numerical 
integration of the (non averaged) equations of motion is required. 



7. Discussion 



In this paper, we provided a framework to describe the bulk properties of the remnant 
of a BHB merger. Our framework is based on PN scaling and fitting the results of 
full numerical simulations. The new formulae are physically motivated, incorporate 
the correct mass ratio dependence, and account for the radiation of angular momentum 
both parallel and perpendicular to the orbital angular momentum. These formulae have 
a symmetric dependence on the mass ratio and spins, while still including the correct 
particle limit. We also extended the successful recoil formula ^ by adding nonleading 
terms that include all the linear dependence in the spins, as well as higher mass ratio 
powers. 

Unlike in the formula for the remnant recoil case, the energy lost by the binary 
during the inspiral phase is a non-trivial fraction of the total radiated energy (and is, in 
fact, the dominant contribution in the small mass ratio limit). We thus included both 
the instantaneous radiative terms in ^ and the binding energy at the ISCO Eq. (g 
(to take into account the secular loss of energy from very large distances down into the 
merger and plunge regime). Similarly, in order to model the final remnant spin, we 
need to take into account both the angular momentum of the system near the ISCO, 
see Eq. ([7]), and the subsequent loss of angular momentum in the final plunge (which is 
particularly important in comparable-mass mergers). 

Using the fitted coefficients in the above formulae, we find that for equal-mass, 
non-spinning binaries, the net energy radiated is 5% of the total mass and the final spin 
is a ~ 0.69, both in good agreement with the most accurate full numerical runs [55] . 
For maximally spinning BHBs with spin aligned and counteraligned we estimate that 
quadratic corrections lead to radiated energies between 10% and 3% respectively. As 
for the magnitude of the remnant spin, the linear estimates are between 0.97 and 0.41 
respectively, with quadratic corrections slightly reducing those values. These results 
show that the cosmic censorship hypothesis is obeyed (i.e. no naked singularities are 
formed) and are in good agreement with earlier estimates [7j. 

The set of formulae ^ and ^ with the fitting constants determined as in the Sec. IS] 
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Table 1. The following parameters give the current best estimates for the constants 
in Eqs. Q and (|6|. These parameters were used to generate the spin-magnitude 
distribution in Fig. [T] 
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can be used to describe the final stage of BHB mergers in tlieoretical, N-body, statistical 
studies in astrophysics and cosmology [MllMlESllMlEZlIMlESlEniEIlESlESllMI^ 
by choosing a distribution of the initial intrinsic physical parameters of the binaries 
{q, 81,82) and mapping them to the final distribution of recoil velocities, spins and 
masses after the mergers. As an example of such an application of the above formulae, 
we calculate the expected distribution of spin magnitudes of astrophysical supermassive 
and intermediate-mass BHs (which are expected to have undergone several mergers). To 
do this, we first consider a set of 10^ binaries with uniform distributions of mass ratio 
(from to 1), uniform orientations of the spin directions, and uniform distributions 
of spin magnitudes. We then use our formulae (see Table [l] for the values of the 
constants that we used) to predict the spin-magnitude distribution of the merger 
remnants and repeat the calculation, again with uniform distributions in mass ratio 
and spin directions, but with this new spin- magnitude distribution (see also [Ml [T2j). 
The resulting spin-distribution, after each subsequent set of mergers, approaches a fixed 
distribution. The spin distribution that results after ten generations of mergers is 
shown in Fig. [TJ The final results are insensitive to the initial distribution and quickly 
converge, in a few generations of mergers, to the displayed curve, which represents a 
universal distribution of the intrinsic spin magnitudes (with a maximum near 0.7 and 
mean in the range (0.5, 0.8)) of the remnant BHs of dry BHB mergers (when neglecting 
accretion). In order to provide a simple analytical model for this distribution, we fit 
it to the Kumaraswamy functional form [67] f{x]a,b) = a62;"~^(l — x")*~^, and find 
a = 5.91 ± 0.04, b = 5.33 ± 0.07. We choose this functional form because it allows for a 
skewed distribution (and fits the results better than a beta distribution), however, the 
fit underestimates the probability for producing small spins. 

We have also considered the effect of wet mergers on the final spin and recoil 
velocities. A first account of accretion effects during the long inspiral phase of 
binary black holes have been given in [68] using the smoothed particle hydrodynamics 
approximation (SPH). To evaluate the accretion effects on the statistical distributions, 
according to [68j, we have considered distributions that at merger have 0.3 < ctj < 0.9 
and orientations within 10 deg and 30 deg for cold and hot accretion disks respectively. 
We have also assumed a flat distribution in mass ratios in the region < g < 1. The 
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Figure 1. The spin magnitude distribution for dry mergers. We plot the distributions 
of spins a = S/m^ of the final remnant after many mergers. This distribution does 
not change significantly following additional mergers and peaks at a « 0.73. We also 
display the distributions representing wet mergers for hot and cold accretion disks. 
They are highly peaked distributions at around a « 0.88 and a « 0.9 respectively. 



results for the final spin distributions are displayed in Fig. [T] and show the dramatic 
change in the spin distributions due to accretion. Note that this accretion effects on 
spin will be less important on black holes with masses larger than lO^M© [5^ . 

The same statistical analysis can be made with the magnitude of the recoil velocity 
of the remnant final black hole when we consider a set of 10^ binaries with uniform 
distributions of mass ratios in the range < g < 1. For dry mergers we consider uniform 
orientations of the spin directions, and uniform distributions of spin magnitudes. We 
evaluate Eq. ^ each time and obtain the distribution with the extended tail beyond 
1000 kms~^ in Fig. [2j The other two distributions correspond to the wet mergers with 
0.3 < ttj < 0.9 and orientations within 10 deg and 30 deg for cold and hot accretion 
disks respectively according to Ref. p8]. The results show a tighter distribution around 
low recoil velocities for cold than for hot accretion disks around the merging black holes. 

Finally, another use of the remnant formulae can be found in modeling waveforms 
in the intermediate and small mass ratio limits using the techniques of Ref. [70] by 
providing an accurate a priori estimation for the background black-hole mass and spin. 
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Figure 2. The the recoil magnitude distribution for dry mergers displaying a tail 
extending beyond 1000 kms^^. We also display the distributions representing recoils 
for wet mergers for hot and cold accretion disks. The cold disk leads to a recoil 
velocity distribution highly peaked at around w sa 80 kms~^ while the hot accretion 
disk extends the magnitude of the recoil to several hundred kms~^. 
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